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Walsh functions form an orthonormal basis set consisting of square waves. Square 
waves make the system well suited for detecting and representing functions with 
discontinuities. Given a uniform distribution of 2? cells on a one-dimensional ele- 
ment, it has been proven that the inner product of the Walsh Root function for 
group p with every polynomial of degree < (p— 1) across the element is identically 
zero. It has also been proven that the magnitude and location of a discontinuous 
jump, as represented by a Heaviside function, are explicitly identified by its Fast 
Walsh Transform (FWT) coefficients. These two proofs enable an algorithm that 
quickly provides a Weighted Least Squares fit to distributions across the element 
that include a discontinuity. The detection of a discontinuity enables analytic rela- 
tions to locally describe its evolution and provide increased accuracy. Time accurate 
examples are provided for advection, Burgers equation, and Riemann problems (di- 
aphragm burst) in closed tubes and de Laval nozzles. New algorithms to detect up 
to two CO and/or C1 discontinuities within a single element are developed for ap- 
plication to the Riemann problem, in which a contact discontinuity and shock wave 
form after the diaphragm bursts. 


I. Introduction 


Accurate simulations of viscous, hypersonic flow generally require high quality grids (hexahedra 
and prisms) aligned with captured shocks and shear layers.! While there is some evidence that 
simplex elements (tetrahedra in three-dimensional space) may be acceptable in the context of high- 
order, discontinuous Galerkin method,? best practice still emphasizes that alignment of cell faces 
parallel or orthogonal to shocks, contact surfaces, and shear layers are less likely to distort local 
flow features with adverse effect of computation of gradient-based quantities (heating and shear). 
Various flavors of multidimensional reconstruction address some limitations of poorly-aligned grids 
in the vicinity of shocks and contact discontinuities, but these formulations are still not as accurate 
as well-aligned, hexahedral and/or prismatic grid systems.*+°® 

Grid alignment becomes more problematic in an unsteady environment. Turbulent boundary 
layers with hypersonic edge conditions are characterized by shocklets that rapidly form, evolve, and 
disappear in the vortical structures that transfer mass, momentum, and energy across the layer.’ 
In like manner, supersonic jets interacting with shock-dominated external flow are characterized by 
rapid, pulsing movement that severely impede the ability to generate prismatic grids that remain 
aligned to all of the interacting discontinuities.* Best practice in these environments is to align 
a prismatic grid or fit a moving boundary to the strongest shocks in the domain and then apply 
well-tuned shock detectors and limiters to weaker shocks that must be captured, regardless of the 
local quality of grid alignment. 
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Is shock capturing with poorly-aligned grids a necessary evil for some, challenging simulations? 
Localized grid refinement is thought to compensate for associated dissipation errors, but is there a 
better, more accurate approach? Algorithms to enable detection of discrete features in a simulation, 
with sufficient accuracy to quantify their evolution across the domain is explored in this paper. In 
this context, it is not enough to simply detect a feature. Sufficient detail regarding the geometry 
of the feature must be extracted so that supplementary equations describing feature evolution may 
be applied. In the case of simple advection, if a discontinuity is detected, it is translated to its new 
location after time At. In the case of a shock, its pre- and post-shock conditions must be uniquely 
defined to enable its movement without smearing using a Riemann solver.? 

In the ideation of an infrastructure for detection and representation of discontinuities, it is 
logical to examine sets of basis functions that are inherently discontinuous. Walsh functions!? 
form such an orthonormal basis set consisting of square waves that are well suited for such fea- 
ture detection and geometry extraction.!* In much the same way that Fast Fourier Transforms 
(FFT) are used to identify frequency distribution, the Fast Walsh Transform (FWT) can be used 
to identify features within the distribution. Indeed, the FWT is commonly used in image processing 
applications.!° Recent papers on the use of Walsh functions in the solution of nonlinear partial 
differential equations (PDEs) have emphasized the ability to derive the solution in transformed 
space.!4>.16 Fundamental operations (multiplication, division, differentiation, and integration) 
are easily programmed. The multiplication of any two Walsh functions g;(x) and g;(x) with indices 
i,j < 2? is another Walsh function g,(a) with index k < 2?. This closure property under mul- 
tiplication enables Walsh function solutions of nonlinear systems of partial differential equations 
(PDE). However, these Walsh function based solution algorithms showed no significant advantage 
over the conventional, upwind-based finite-volume techniques they were intended to replace.® The 
important lesson from these tests is that, while Walsh functions are well suited to representing 
discontinuous functions, the discontinuities may still be dissipated if the discretization algorithm is 
not formulated in a way that preserves the detected features. 

The vision for the FWT fit is to work like a shock-capturing method with respect to detection 
but work like a shock-fitting method with respect to use of more accurate analytic relations. The 
long term goal of this work is to enable accurate simulations of viscous, hypersonic flow on simplex 
elements (triangles in 2D, tetrahedra in 3D). For steady problems, the objective is accuracy that 
meets or exceeds that obtained on a well-aligned prismatic grid system. For unsteady problems, 
the objective is to create an infrastructure that accurately detects rapidly evolving shocklets and 
computes their movement across the domain using analytic jump conditions. The goal requires: 
(1) the ability to compute viscous quantities (heating, shear) on simplex elements that are not 
generally well-aligned with the shear layer, and (2) the ability to detect shocks and maintain 
proper jump conditions on poorly-aligned grids. Much progress has been realized toward the first 
goal through a hyberbolic formulation of the governing equations that enables highly accurate 
resolution of gradients on simplex elements.'’ Accuracy has been demonstrated on such grids 
with randomized perturbations that intentionally disrupt grid quality.!® This paper is the second 
of two to address the detection and tracking of discontinuities with Walsh functions. The first 
paper!® developed the mathematical foundation for using Fast Walsh Transforms (FWT) to detect 
discontinuities while simultaneously providing a baseline polynomial fit to smooth regions spanning 
a computational element. The most fundamental aspects of feature detection are developed there 
for one-dimensional simulations. This paper reviews highlights of the first and then presents new 
material on the time accurate detection and tracking of discontinuities for the Euler equations. 
It essentially seeks to replace shock capturing with shock fitting.?°?! Preliminary investigation 
suggests that multidimensional extensions could utilize a multidimensional Walsh function series 
using the same fractal recursion formulation on simplex elements that is already defined for the 
one-dimensional series. 
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This paper is organized as follows. Walsh function fundamentals are reviewed in I. Highlights 
from the first paper!® describing feature detection and algorithm formulation using Walsh functions 
are reviewed in Sec. II. The theorem addressing Walsh Root Function annihilation of low-order 
polynomials is discussed here. It becomes the cornerstone for a Weighted Least Squares Fit based 
on the FWT, with results presented in Sec. IV. Feature detection in the presence of several off- 
nominal conditions, such as noisy distributions or two discontinuities, is discussed in Sec. V. 
The principal focus of this paper is to enable preservation and tracking of discontinuities that 
develop in the solution of a time-accurate PDE by exploiting the FWT capabilities for feature 
detection. To this end, a conservative formulation of one-dimensional equations are presented in 
Sec. VI. The key innovation here is to formulate the flux relative to the FWT curve fit, including 
the discontinuity. All corrections to the baseline fit are then based on a Cp continuous baseline. 
The new algorithm is tested through simulations of time-accurate advection, Burgers’ equation, 
and quasi-one-dimensional nozzle flow and shock tube flow in Sec. VII. Concluding remarks and 
directions for future research are finally presented in Sec. VIII. 


II. Walsh Function Fundamentals 


Let gn(x) be a basis function over the domain xq < x < 2» with n segments. The orthonormal 
basis function requirements” that 


- On Eon dE = Onn (1) 


enable a recursive algorithm to define gn(az) based on the following constraints. Constrain |gn(x)| 
equal to a constant across the entire domain (i.e., only allow square waves of constant amplitude). 
Let index p > 0 define a segment size dzp. 


dry = (ay — ta)/2 (2) 


An additional constraint on segment lengths within g,,(x) requires that at most two segment lengths 
(dxp and 2dz,) are allowed to create n segments spanning the domain. These new constraints allow 
basis functions to be grouped according to the smallest segment size. A basis function gn (x) belongs 
to group p if 
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Figure 1. Basis function in group p=0. 


As an example, consider the evolution of the first four groups of g,(x) corresponding to0 <p <3 
on the domain with x, = 0 and 2, = 1. For p = 0, the only element of the group according to Eq. 3 
isn = 1. A single segment with length given by Eq. 2 spans the domain. Therefore, gi(a) = 1 (Fig. 
1) is the first basis function, belonging to group p = 0 and satisfying Eq. 1. (While gi(a) = —1 
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Figure 2. Basis function in group p=1. 


also satisfies Eq. 1, it is convenient to adopt a convention where the sign of the function on the 
first segment is positive.) 

For p = 1, the only element of the group according to Eq. 3 is n = 2. Two segments with 
length equal to 1/2 span the domain. Therefore, 


1 for O<%¢22 
g2(x) = ‘ ; 
—1 for Si | 


where function go(x) (Fig. 2) is the second basis function, belonging to group p = 1 and satisfying 
Eq. 1. The value of the function at interior segment boundaries is set to 0, the average of the 
function on either side of the boundary. (The vertical line through interior segment boundaries 
in Fig. 2 provides strong visual impact of the discontinuity and is not meant to imply all values 
between —1 and 1 are included.) 

The group defined by p = 2 is the first to have more than one function, g3(2) and g4(x). The 
allowed segment sizes in this group are 1/4 and 1/2. Each function is initialized with the smallest 
segment at the beginning and end of the domain. For n = 38, this initialization leaves a single 
segment of length 1/2 available to span the interior of the domain. For n = 4, the initialization 
leaves two segments of length 1/4 to span the interior of the domain. The corresponding functions 
are shown in Fig. 3. By inspection, it is clear that these first four basis functions satisfy Eq. 1. 


1 for O<x<} 
1 for O<x<} : 
7 " —1 for f<u<s 
g3(@) = 4-1. for ee a eer ga(z) = : 
1 for 5<a<# 
1 for #<a<1 
-1 for 3<a2<1 


The group defined by p = 3 includes four functions, g5(a) through gg(x). The allowed segment 
sizes in this group are 1/8 and 1/4. Each function is again initialized with the smallest segment at 
the beginning and end of the domain. This is the first group where the distribution of segments 
required to satisfy Eq. 1 for all combinations of functions in groups 0 — 3 becomes nontrivial. 
Still, the number of permutations of segment size distribution is small and a solution, evident by 
inspection, is shown in Fig. 4. 

A fractal-based recursion algorithm is used to define the wave pattern for ever larger numbers 
of segments.!® A function f(x) may now be represented by an infinite series 


Ss" Angn (x) (4) 
n=1 
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(a) g3(x) (b) ga(x) 


Figure 3. Basis functions in group p=2. 
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(c) gr(«) (d) gs(x) 


Figure 4. Basis functions in group p=3. 


over the interval 7g < x < x». The coefficients A, are computed 


An = ‘il FG, ade (5) 


Note that, for domains comprised of a maximum N = 2? segments, the Fast Walsh Transform 
(FWT)!® will be used to compute the wave components A, from f, in order pN operations. If f, = 
{eke (x)dx (integrated average), then a continuous representation is defined. If f, = f(an), 
then a discrete representation is defined. In the derivations that follow, a discrete representation 


is assumed. 


III. Feature Detection from the Fractal Fingerprint 


The fractal nature of the Walsh basis functions!® introduces repeated patterns in the wave 
component magnitudes that can be exploited for many applications. These patterns are accentuated 
in a mapping defined as the fractal fingerprint (FFP). The FFP has been presented in detail 


5 of 29 


American Institute of Aeronautics and Astronautics 


previously! and highlights from that work are summarized here. 


Walsh functions in group p are simply step functions with a baseline step size equal to 
(ap — Zq)/2? and a maximum of two adjacent steps with the same sign. A mapping that 
retains group identity for the plot of A, in the FFP is defined: 

m=2?+1—n for 2 l<n<2? 


(6) 


For singularity-free functions, a recursive algorithm is evident that enables prolongation - an 
estimation of all but one element, r(p), of a group p coefficient from a given group p—1. The 
new, unique element r is defined as the root index and the Walsh function with root index r 
is defined as the Walsh root function g,(p)(2). 


It is proven!® that the inner product of the Walsh root function Jr(p) (x) with any polynomial 
of degree p—1 over the interval xg < © < Zp is identically zero. This proof enables an efficient 
algorithm for a least squares polynomial fit through data sets spanning 2? cells in the interval. 


Coefficients of the polynomial curve fit, 


f(a) = MS, apn 
k=0 


to the distribution f(x;) for 0 <7 < 2? are computed from 


(7) 


B,0),0 Bro) Bro),2 B,o),m ao A,(0) 
0 Baya Bray. Bray aa ay A,(1) 
0 0 Bre),2 B,(2),m a2 | =| Ap) (8) 
0 0 0 By(m),m am Ae (ea) 


where A,(;,) are the coefficients of the Walsh root functions obtained through the Fast Walsh 
Transform (FWT) of f(x;) and B,(,)m are constants, equal to the coefficients of the Walsh 
root functions obtained through the FWT of 2?’. 


Given a distribution, f(x;), and a shifted distribution f(a; + Az), across 2? cells on an 
interval x, < x < 2», two essential properties of their FFPs can be identified to enable 
feature detection of a discontinuity with a jump magnitude greater than |h|. 

— The average of |.A,,| for 2?-' < n < 2? must be greater than |h|,/xp — %q/2? 

— The pattern of sign changes of A, for 2?-! < n < 2P determines the location of the 


discontinuity within the interval 


A distribution with a Cp discontinuity may be approximated by 


f(x) = So one” +hH (a, x*) 


k=0 


(9) 


where 
0 for 


1 for 


< * 
a (10) 
cr > x 


The coefficients a, are determined by subtracting the influence of a detected shock from a 
baseline shocked distribution. 
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— First, compute and store the constant Walsh Root function coefficients Brow) , for 


H(a;,x*) for every ordered pair (xz, 2,41) bounding a shock location «* occurring in 
the interval (xq, Xp). 


QP 
Bris a: = Soe Ae tee Kt ee Coarse 11) 
i=l 


— Second, the Walsh Root coefficients for a smooth distribution A,.,,) are now computed 
in Eq. 12 from the Walsh Root coefficients of the shocked distribution Arc) minus the 
rescaled Walsh Root Function coefficients for the Heaviside function Brow) , with the 
shock location given by index k. 

Ar(y = Ar(ny — RB r(x) k (12) 


r(K 
IV. Polynomial Curve Fits from the Walsh Transform 


Representative curve fits for several functions are presented in Fig. 5. The discrete data f; 
being fit are shown as black squares. Curve fits through the data are shown as a black dashed 
line for the 4th degree polynomial (m = 4) and as black solid line for the 6th degree polynomial 
(m = 6). Both lines are contained within the symbols for each example. The difference between 
the discrete data and the curve fits are shown on the right scale with red lines. The difference is 
identically zero for f; = x? in Fig. 5a because the algorithm in Eqs. 8 will exactly recover the 
polynomial coefficients if f; is formed from a polynomial of degree less than or equal to m, the 
degree of f(x;). Errors associated with lack of precision begin to manifest at m > 8 - the Walsh 
Root Function coefficient associated with r(9) = 342 tends to have a very small magnitude (related 
to the sum 2° differences of nearly equal quantities) and this loss of precision directly affects the 
calculation of the leading coefficient ag. 


0.0001 


'-5E-05 


(a) fi=ai (b) fi = cos(2ras) (c) fi = exp(ai) 


Figure 5. Comparison of polynomial curve fit f(x;) = rip ex} (black lines) to singularity-free func- 
tions, f; (square symbols) with 1 < i < 64 (p = 6). Their difference is shown on the right scale for 
polynomials of degree 4 (dashed red line) and degree 6 (solid red line). 


The value of fj41/2 at the interface of two cells will be required in subsequent simulations 
in Sec. VI. To assess accuracy of interface interpolation, formulations of the interface value are 
compared to the special case where an exact interface value f.;+1/2 is available. Cell interface 


values will be computed on the basis of the fit value at the interface, f (2;41/2), and interpolated 
corrections from neighbor cells, (fi — _ (a;)). Consequently, the interface value is a function of the 
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curve fit polynomial f (x) of degree m that is fit to 2? discrete points where m < p. A third-order 
interpolation with this property is defined below in Eq. 18. 


fisya = F (titi) + : 3 (fin = F(ais1)) eo (F - F(xi)) = (fi-2 “ F(ai-2))| 
= ; (3fia + 6f; = fia) = ; (Fait) 7 2f (@i41/2) + F(xi)) 
+ & (Flea) —2fe) + Fle) ard order (13) 


The error norms (£1) of the corrected interface values from Eq. 13 are defined in Eq. 14. The 
error norms of a reference, uncorrected interpolation formulation (L1,e¢) are defined in Eq. 15. 


2P—1 
Iy = > |fittye— feit1a| Av (14) 
=i 
2P—1 
Ly ref = Ss" 8 (Sie OR — fA) —Fegiajo| Ae (15) 
=] 


Note that 2441/2 = (wit+ai+1)/2 and fei41/2 is the interface value for one of three test functions 
(exp(aj41/2), Cos(27xj441/2) or 1/(1 + 2441/2). The last function here replaces the monomial ies 
because the Walsh transform fit identically reproduces this function. 


slope = -3 
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(a) fi = exp(xi) (b) fi = cos(27r2;) c) fe =1/(1 +2) 


Figure 6. Variation of £; error norm with p for 3rd-order accurate interface interpolation on the 
interval 0 < x < 0.5 with 2” discrete points comparing conventional interpolation (black line) to Walsh 
Transform curve fit corrected interpolation (red lines) for polynomial fits varying from degree 3 to 
degree p (3 << m< p). 


It is instructive to compare the error norms of the corrected interpolation formulations from 
Eq. 14 to the error norms of the reference, uncorrected interpolation formulation in Eq. 15 in Fig. 
6. In these figures, each increment in p corresponds to a factor of two decrease in the mesh size. 
A slope of -3 in Fig. 6 indicates that the error norm drops by a factor of eight for an increment 
in p equal to one. All of the interpolation formulations are third-order accurate except for the 
m > 6 in Fig. 6a in which the formulation bumps up against machine zero capability with double 
precision arithmetic. Orders of magnitude increase in accuracy are demonstrated as the order of 
the polynomial fit m increases for the correction terms. Note that in Fig. 6b the variation of error 
norms is nearly identical for m = 3 or 4, m = 5 or 6, and m = 7 or 8. This behavior derives from 
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the fact that the even powers of x in the fit to cos(27x) are almost identically zero on the interval 
0 < x < 0.5. Consequently, the next even power of the curve fit provides no new information 
beyond the previous odd power. 

The rate of improvement with increasing m depends on the quality of the curve fit. The red 
lines are generally observed to lie below the reference black lines for smooth, monotone functions. 
Even if the function has a local maximum or minimum, improvement is usually noted unless an 
attempt is made to capture distributions with multiple extrema. 

Finally, the operation count required to compute the correction terms is a relatively small 
overhead with respect to that required for the reference interface interpolation. Evaluation of the 
polynomial coefficients in Eq. 8 requires order p2? operations; however, this single fit is used to 
compute the correction to 2? — 1 interface values. 
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(a) original distribution and curve fit, h = 1.5 
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Figure 7. Feature detection within the distribution given by Eq. 16 with discontinuity at 111 < x2* < x12. 


The problem of detecting a shock in the context of a more general, underlying distribution is 
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exemplified in Fig. 7 with distribution defined in Eq. 16. 


1 


fi, = cos(27a;) + hH(2;, x*) for U<a;<1, 2 == -)Agv, Ag= 6A 


(16) 
Examples of curve fits f(x;) including the detected shock location are provided with solid black 


lines in Figs. 7a and 7c. The corresponding FFPs used to derive these fits are shown in Figs. 7b 
and 7d. 


V. Special Cases 


A. Noisy Distributions 


2.5 ; 
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(a) noisy distribution and curve fit with m = 3 (b) FFP of noisy distribution 


Figure 8. Curve fit and FFP of distribution given by f; = cos(27z;) + 0.2(rnd; —0.5)2 on interval (0,0.5). 


Random numbers with magnitude 0 < rnd; < 1 are mapped to the function f; = cos(272;) 
on the interval (0,0.5) in order to test the algorithms for determining polynomial fits and feature 
detection on noisy distributions. If m > p—1 in Eq. 7 and the distribution f; is noisy, then it 
is usually observed that the resulting curve fit using the algorithm defined by Eq. 8 is poor. The 
magnitude of the highest-order polynomial coefficient a, in such a case is usually too large and the 
fit oscillates through the distribution with large excursions above and below the mean. It is found 
that m < 4 through distributions with p > 6 produce good curve fits through noisy data, as shown 
in Fig. 8. The noise is effectively filtered by basing the curve fit on the truncated Walsh function 
transform. 

The ability to detect a shock in the presence of a noisy distribution is demonstrated in Fig. 
9. The two randomly generated test cases presented in Fig. 9 highlight the challenge of detecting 
a shock in the presence of noise with amplitude equal to 0.2. In most cases, a result consistent 
with Random Sample 1 in Fig. 9a is obtained. The pattern match in the FFP (Fig. 9b) is 
somewhat disrupted by noise; however, the disruption was not sufficient to thwart the detection 
algorithm. Occasionally, the randomly generated data produces an FFP with too many disruptions 
to an accepted pattern to identify a shock location, as shown in Fig. 9c and d. A similar test in 
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NO 
OAvd bom 
Tc 
" 

On SBwWNn= 8 


0 0.2 0.4 0.6 0 8 16 24 32 
x m(n,p) 
(c) Random Sample 2: noisy distribution with shock (d) Random Sample 2: FFP of noisy distribution with 
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Figure 9. Curve fit and FFP of distribution given by f/f; = cos(27x;) + H(ai,x*) + 0.2(rnd; — 0.5)2 on 
interval (0,0.5) with 211 < x2* < 212. 


which the amplitude of the noise was reduced to 0.1 encountered no failures of the shock detection 
algorithm. 
B. Distributions with Two CO or C1 Discontinuities 


Feature detection in a hypersonic flow solver must be able to accommodate the interaction of two 
shocks in the same element. An analysis of the FFP for two Heaviside functions reveals that 
the resulting distribution of Walsh coefficients is nothing more than a linear combination of each 
individual function’s coefficients. A bilevel distribution of the group p coefficients is observed, from 
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which it is possible to reconstruct the effect of each Heaviside function individually. A representative 
curve fit is shown in Fig. 10a. 

In like manner, the detection of C1 discontinuities is accommodated by taking the derivative of 
the transform!® and then search as before for the CO discontinuities. The differentiation operator 
is simply a matrix multiply of the Walsh coefficients. An example with C1 (slope) discontinuities 
is generated with H'(x,x*) = (a — 2*)H(x,2*) and is presented in Fig. 10b. 
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0 0.2 0.4 0.6 
» 4 
(a) fi = cos(27x;) + 0.8A (x, xj) — 0.5H (2, 73) (b) fi = cos(2rax;) - 20H'(x, xt) + 20H'(x, x3) 


Figure 10. An example of feature detection for a two discontinuity distribution with x29 < xj} < x30 
and 2733 < © < #34 for p=6 and m=4. 


VI. Formulation Relative to the FWT Fit 


The FWT curve fit with feature detection provides new opportunities for the solution of differ- 
ential equations. Simple examples are considered in this section. The essence of the approach is to 
linearize the system of equations with respect to the FWT fit. 

The value of any variable calculated from the reference FWT fit is indicated with a tilde over the 
variable name. The formulation of any interface value at index i + 1/2 between indices 7 and i+ 1 
requires a single curve fit that spans all cells with an explicit influence on the interface value. Each 
FWT is evaluated at the beginning of a timestep. A very simple way to satisfy this requirement 
in one space dimension is to run an FWT fit centered at each cell interface. Such an approach 
is not readily extended into multidimensional simulations using simplex elements. In the present 
formulation, a domain is divided into elements. (See Fig. 11.) Each element contains N = 2? cells. 
A single FWT curve fit is derived for each dependent variable spanning the element. In general, 
a fourth-order accurate formulation of inviscid flux at cell face i + 1/2 requires support from cells 
i —2toi+3. The interior cell walls from N/4+1 to 3N/4, indicated by the upper black brackets 
in Fig. 11, are all defined relative to the elemental FWT fit. The remaining cell walls within N/4 
cells of an elemental boundary, indicated by the lower red bracket in Fig. 11, are defined relative 
to a blended FWT fit spanning the interior elemental boundaries. 

The model equation system is given by 


Oq Of Oh = 
mt act dy 1 8aet) =9 (17) 
12 of 29 


American Institute of Aeronautics and Astronautics 


Element 1 Element 2 
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Figure 11. Schematic of the stencils used to define the inviscid flux across cell boundaries. Each 
element is composed of 16 cells. An FWT fit is applied across each element. The upper black bracket 
defines the cell boundaries derived from the elemental FWT fit. A blending FWT fit is applied across 
element boundaries. The lower red bracket defines the cell boundaries derived from the blended FWT 
fit. 


A fourth-order Runge Kutta method is used to advance the solution in time.?? Let 


F(x, t, q) = Ox a Ox? 


Then, given an initial distribution at time t = to, the solution at t = to + At is given by 


+s(q,x,t)|. (18) 


a(x, to + At) = ale to) + 5 (M(x) + 2M2(2) + 205(0) + O4(2)) (19) 
where 
HG) = AW e@h,a@%)) (20) 
Q(z) = AtF(2, to + sat, a(c, to) " 5% (2)) (21) 
Q3(x) = AtF(2, to + sat. alc, to) + 5(2)) (22) 
WG) = Aig Abatedys- ase) (23) 


A conventional, fourth-order accurate formulation in space for the inviscid (or viscous) flux 
terms are defined by 


(5) = (27 (fi41/2 — fi-1/2) — fi43/2 — f;_3/2)) /(24Az) (24) 


When discontinuities are detected in the FWT they cannot be located with any greater precision 
than one cell size. The detection algorithm simply informs that a discontinuity has been detected 
between cell center 7 and cell center i+ 1. For convenience, the discontinuity is assigned a location 
at the cell interface, i+1/2. In this case, the limiting value of f when approaching the cell interface 
from below, Fiat Jo is not equal to the limiting value of f when approaching the cell interface from 
above, ue jo" It also follows in this case that the variation of f between 4 /2 and «71, 2 is 
continuous. Equation 24 is therefore modified to account for the possibility of a discontinuity at a 
cell interface as follows. 


at 7 7 - 
(5) = (26 Oye Baye) — Cae Bhape) — Cae Btgn)) (24dz) (28) 
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Even though fourth-order spatial accuracy is lost in the presence of a discontinuous jump (except 
for some special cases), it is expected that this formulation will preclude the need for additional 
limiting or adjustment of stencils. This formulation is used for the advection and Burgers’ equation 
test cases that follow. A second-order accurate formulation is expressed by 


&) = (faa hae) M2) a 


will be defined relative to the FWT fits in Sec B. Viscous flux terms 


Inviscid flux terms fi /2 
hj+1/2 are defined in C. 


A. Preservation of Discontinuities 


The right and left state of a discontinuity in an element may be defined in one of two ways. If the 
FWT is a good fit to the data spanning an element, then the right and left states across a cell wall 
z* = 2441/2 are computed from Eq. 9 with x = 2* + eAz, where € = 0.0001 herein. A good fit is 


defined as - z 
von=i 19? (vi) — 9? (ai) 

Iy= * u - . | < Ly target (27) 
where Lj target = 0.01 for the Euler equations. Note that the Euler equations use FWT fits of the 
primitive variables, q?. (This conditional was never needed for the simpler advection and Burgers’ 
sample problems.) If the FWT fit fails this test, then an exact Riemann solver? is applied across 
every cell wall in the element and the right and left states are defined by the values at adjacent cell 


centers. 


expansion contact 


fan \ \ shock , 


p,u 


Figure 12. Riemann problem schematic example: given an initial state 1 and state 4 on either side 
of a diaphragm, three waves are formed upon the abrupt removal of the diaphragm. These waves 
include a shock or expansion fan moving to the left separating intermediate state 2 from state 1 and 
a shock or expansion fan moving to the right separating intermediate state 3 from state 4. A contact 
discontinuity separates intermediate states 2 and 3. The variation of pressure, density, and velocity 
are shown, with pressure and velocity continuous across a contact discontinuity. 


In general, a discontinuity in the Euler equations decomposes into three waves as defined by 
the Riemann problem illustrated in Fig. 12. The algorithm for detecting and tracking these 
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discontinuities is represented in the flow chart of Fig 13 with additional details as follows. An exact 
Riemann solver? is used to determine the intermediate states 2 and 3 given the left and right states, 
1 and 4, respectively, from Fig. 12. The Riemann solver will also provide the wave velocities of any 
shock or expansion fan separating states 1 from 2 (left running) and 3 from 4 (right running) as well 
as the contact surface velocity separating states 2 and 3. For the purpose of feature detection in 
this work, a detection threshold is defined for a shock as a function of the pressure ratio across the 
wave, (Pnigh/Plow)thr- In like manner, a detection threshold is defined for contact discontinuities as 
a function of the density ratio, (ppign/Plow)thr- Expansion fans are not preserved as discontinuities. 
Feature tracking requires that the bounding states and velocity of the leftmost detected wave be 
stored as q;. qi, inner’ 2nd Vy. In like manner, the bounding states and velocity of the rightmost 
detected wave are stored as dp: GR inner’ and Vp. A key element of this algorithm is to accept 
the fact that, even though a wave grouping as represented in Fig. 12 may be comprised of distinct 
waves after some small amount of time, it is not possible to observe that distinction in discretized 
space until at least one element of the group moves past a cell center. Consequently, an FWT 
detected discontinuity between any two cell centers must include the leftmost and rightmost states 
given by qi, ee VL, dh, nian and Vr. The detected discontinuities identified as a shock 


Evaluate new cell center value 
after discontinuity passes over 


Execute the FWT fit. 


(1) Evaluate exact Riemann solver 
across FWT detected discontinuities. 
(2) Preserve discontinuities that meet 
detection threshold. 

(3) Use Roe flux approximation across 
discontinuities below the detection 
threshold. 

(4) Use linearization wrt to FWT 
across all other interfaces. 


L 
(1) Evaluate exact Riemann (1) Advance solution at all Any 
solver across all interfaces. cell centers one time step. ARVN IAN 
(2) Preserve discontinuities (2) Calculate movement of past cell 
that meet detection threshold. all discontinuities. aan 
(3) Use Roe flux approximation — 
across all other interfaces. 


Figure 13. A flow chart for the detection and tracking of discontinuities. The logic proceeds from 
the upper left corner and is executed on every element in the domain. One cycle through the chart 
corresponds to one time step. Options in red text are only required for the Euler equations. 


or contact surface serve as internal boundary conditions for the flow solver. Therefore, the finite 
volume formulation for any cell bordering a discontinuity will use 


fi = f(qz) (28) 
fi 4/2 = f(qp). (29) 


The detection algorithm starts afresh every time step. Shocks will appear or disappear based 
solely on the detection parameters!? used within the FWT and the threshold parameters, 
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(Phigh/Plow)thr 20d (Prigh/Plow)thr, filtering the Riemann solver. Still, it is necessary to maintain 
information regarding the previous location of a detected discontinuity so that its movement across 
the element may be computed. To enable shock movement, a data storage entity defined as a 
floater is created. Two floaters per element are allowed because the FWT (as currently configured) 
can only discern at most two discontinuities. Floaters are initialized to a null state. (See Fig. 14.) 
A newly detected discontinuity is assigned to a floater to enable its location s relative to 2j41/9 
to be tracked. A detected discontinuity matches an existing floater location if it lies between the 
same two adjacent cell centers. If, after a new time step, a floater is not matched, then the feature 
disappears as a discontinuity and the cell interface is no longer treated as an internal boundary 
condition. Though not shown in the flow chart, floaters are free to cross element boundaries. 


=O If s, < -Ax/2 or sp > Ax/2 Update the discontinuity 
Initialize floaters to then move the floater ___ location within the floater: 
null state and apply one cell to left or right, sMi= 5" +V, At 
FWT to element. respectively and update saMtt= sal + Vp At 
— cell center solution. f 
Loop over ie | t=t+At Reset any unmatched 
detected by the FWT Apply FWT to element. floaters to null state. 


(A shock “dies”.) 


Does the shock or 
contact position 

match an existing 
floater location? 


Update Riemann states 
and velocities in the 
matched floater. 


Initialize location to s =0 
and store Riemann states 
and velocities in the floater. 
(A “shock” is born.) 


Figure 14. A flow chart for use of floaters, storage entities within an element that maintain data 
required for shock tracking. Logical flow starts from the upper left for each element. 


In the case of the Euler equations, the FWT may detect a discontinuity that is too weak to 
register as either a shock or a contact surface. In this case, no effort is made to preserve the 
discontinuity as a distinct entity and a Roe’s approximate Riemann solver is used to provide a 
definition of f;41/2. Consequently, the flux for this transitional state between a clearly continuous 
and clearly discontinuous profile assigns, i) j2= fit /2 


B.  Inviscid Flux 


The previous section defines fa 2 for cases where a discontinuity is detected by the FWT between 
cells 7 and i+ 1. In this section, a continuous variation is expected across the interface and the 


FWT is utilized to provide enhanced accuracy as discussed previously in Sec. IV and illustrated in 


Figs. 5 and 6. 
The flux Jacobian is diagonalized in the usual way in Eq. 30. 
Of 
=~ =A=R''AR 30 
= (30) 
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The conversion from primitive to characteristic variables is defined relative to the FWT fit at cell 
centers 7 using the matrix defined in Eq. 31. Differential contributions to the flux emanating from 
left and right running waves are evaluated with respect to the FWT fit at cell interfaces i + 1/2 
using the matrices defined in Eqs. 32 and 33. 


R; = R(G(x)) (31) 
Aisi2 = A(Q(%i41/2)) (32) 
Raye = R(t) (33) 


The difference between the actual value of the dependent variable q; and the FWT fit to that same 
variable at cell center x; is given by Eq. 34. This differential represents a correction to the FWT 
value of the dependent variable at cell centers. The corresponding differential of the characteristic 
variable is defined in Eq. 35. 


ogi = Gi — 4(zi) (34) 

oq ; = R,6q; (35) 
Note that the variation of these differential quantities is expected to be continuous and differen- 
tiable, even across discontinuities. 


A fourth-order accurate interpolation formula is used to bring the cell-centered, differential 
corrections from the left (Eq. 36) and from the right (Eq. 37) to the cell interface at i + 1/2, 


bd/ti1/2 = L(6q'i 41; 64';, 0q/;_1, 0q';_2) (36) 
6d 41/2 = L(6q j, 9G 41,99 142,94 i43) (37) 


where the interpolation formula for equally spaced cells is given by Eq. 38. 
L(a,b, c,d) = (5a + 15b — 5c + d)/16 (38) 


Discerning right and left running corrections to the flux at interface i+ 1/2 employs the eigenvalue 
and eigenvector matrices in Eqs. 39 and 40. The right and left corrections are finally applied to 
the FWT reference of flux f(q(aj+1/2)) in Eq. 41. 


1 


bf) 1/2 = Rist [Aisi/2 oe |Ais1al] bd’ i1/2 (39) 
= ditt a 

Of 170 am a Rist/2 [Ai+1/2 = |Ai+1all 6d i41/2 (40) 

finvyo = £(G(i41/2)) + bf 2 + OFF 172 (41) 


C. Viscous Flux 


Whereas the inviscid flux terms were defined relative to the FWT fit, the viscous terms are defined 
relative to q;. The inviscid formulation is designed to detect and preserve discontinuities. The 
viscous formulation, consistent with the physical nature of dissipation, must retain the ability to 
dissipate (spread out) the shock. The competing influence of these formulations on shock structure 
will be discussed further in the Burgers’ equation examples that follow in Sec. B. The fourth-order 
accurate formulation of viscous flux on a cell wall is therefore defined by 


O 
biyre=¥ (58) | =v(2raiea ~ a8) - (ain2 ~ a2) /24Ac. (42) 
U/ 4541/2 
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This formulation utilizes the identical stencil defined for the inviscid flux in Fig. 11. A modified 
version of Eq. 42 recognizes that the available resolution may be insufficient to resolve a shear 
layer. Equation 43 prevents formation of non-physical extrema across the discontinuity if the 
physical viscosity coefficient is too small. 


biyiy2 = Yo] >, kanee cs + (27(Sai+1 — bas) — (Sai42 — Sqi—1)) /24Ax (43) 
k=1 
+ (v—W)[ (27( git — Gi) — (Gita — ai-1 )) /24Az] 


where v9 = min(v, Ax/2) and the ax are the FWT derived curve fit coefficients of Eq. 9. 


VII. Sample applications 


A. Advection Equation 


The Advection equation is defined from the model system (Eq. 17) with f = cq, h = 0 (v = 0), 
and s = 0 with q a scalar. The domain is given by 0 < x < 1 and a periodic boundary condition 
is applied such that q(0,t) = q(1,t). The wave speed c = 1 is a constant and the matrices in Eqs. 
31 - 33 are scalars, identically equal to 1. Explicit time advancement is used with At = Az. The 
exact solution is expressed with a transformation from x and t to z by 


‘iss x—ct for x>ct (44) 
1 —(\a —ct| —trun(|x—ct|)) for «2 <ct 


so that q(z) = q(z,0). 

Results for the advection of a Gaussian profile are shown in Fig. 15. There are no discontinuities 
to be detected in this example. The computed profile is indistinguishable from the exact solution 
after 20 cycles across the domain in Fig. 15a for a simulation with 8 elements, 2° = 64 cells 
per element (p = 6), using an FWT fit polynomial of degree 6 (m = 6). The Ly error norm for 
several simulations relative to the exact solution after 20 cycles is presented in Fig. 15b, where 
Ly = 95 Aa|q(a;, 0)—q(z;)|. The total number of cells in a simulation equals the number of elements 
times the number of cells per element (2?). A four element simulation with 2? cells per element 
has approximately the same error norm as an eight element simulation with 2?~! cells per element. 
The slope of log,)(L1) with respect to p for simulations with the same number of elements equals 
-4, indicating a fourth-order accurate simulation in space and time. The A symbols in Fig. 15b 
indicate results for which the FWT fit was circumvented by setting q; to zero. This result shows a 
factor two reduction in the error when using the FWT fit (square symbols) as compared to a more 
conventional formulation (delta symbols), even for a case in which there are no discontinuities. 

Results for the advection of a triangular profile that is Co continuous, but has Cj discontinuities 
at z = 0.25, 0.50, and 0.75, are shown in Fig. 16. The computed profile shows slight rounding 
of the corners after 20 cycles in Fig. 16a for a simulation with 8 elements, 2© = 64 cells per 
element (p = 6), using an FWT fit polynomial of degree 6 (m = 6). The corresponding Ly; error 
norm for several simulations relative to the exact solution after 20 cycles is presented in Fig. 16b. 
The presence of C; discontinuities disrupts the fourth-order convergence observed for the previous 
case. The approximate order of convergence is 1.5. Simulations with and without the FWT fit are 
essentially identical - suggesting that discontinuities in slope may be an additional feature worth 
detecting in future work. 

Results for the advection of a cosine profile superposed on a Heaviside function are shown in 
Fig. 17. This is the first example with a Co discontinuity included in the initial profile. The 
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Figure 15. Advection of Gaussian profile f; = exp(—10(z; — 0.5)) after 20 cycles across the domain. 
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Figure 16. Advection of triangle profile f; = 1— 4|z;—0.5)| for 0.25 < z; < 0.75 and f; = 0 for z; < 0.25 or 
z, > 0.75 after 20 cycles across the domain. 


discontinuity should advect to the right with speed equal to one. If the discontinuity is dissipated 
by the numerics in the solution of the advection equation, then there is no restorative action to 
resharpen the discontinuity - all characteristics have wave speed equal to one. Consequently, feature 
detection is exploited here to preserve the magnitude of the detected discontinuity and to move it 
with the wave speed, as outlined in Fig. 14. In this simple advection case, there is only a right 
running wave with wavespeed Vg = 1. When the discontinuity crosses a cell center x;, the cell 
is updated according to the translated FWT curve fit given by Eq. 9. (Note that this update 
algorithm is simpler than that described earlier;'® however, the results are nearly identical.) 

This tracking algorithm effectively preserves the discontinuity through 20 cycles across the 
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Figure 17. Advection of cosine profile with a jump, f; = cos(27z;) + H(z,0.5) for 0< 2 <1. 


domain, as seen in Fig. 17a and b. If the tracking algorithm is not engaged and a conventional 
limiter is not enforced, then ringing across the discontinuity is evident in Fig. 17c. 

The error norms for several simulations varying the number of elements spanning the domain 
and the order of the FWT fit in each element are presented in Fig. 17d. Engaging the discontinuity 
tracking algorithm provides orders of magnitude more accurate results than without tracking and 
provides a convergence rate approaching third-order in space and time. It is important to note 
that the tracking algorithm diverges for m > 5. It appears that the extra degree of freedom in the 
FWT fit enables waviness across the detected discontinuity that amplifies with increasing time. 


B. Burgers’ Equation 

Burgers’ equation is defined from the model system (Eq. 17) with f = 59”, h = pS, ands = 0 
with q ascalar. The domain is given by —1 < x < 1 and fixed boundary conditions are applied such 
that q(—1,t) = 1 and q(1,t) = —1. The initial condition is given by q(x,0) = —zx. The solution 
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evolves to a steady state given by q(x) = C tanh( $2). The constant C' is defined to recover the 


boundary values. In the limit v > 0, q(x) = 1 — 2H(x,0), where the Heaviside function was 
defined in Eq. 10. The matrices in Eqs. 31 - 33 are: R; = [1], Aj+i/2 = [Gi41/2], and R; 


=U), 
i+1/2 [ 
Achieving high-order spatial accuracy in these simulations is simplified by using the value of the 


steady-state analytic solution for any formulation that requires cell center information beyond the 
boundaries. 


In this example, a detected discontinuity between cells i and 71+ 1 has a single wave speed 
approximated by V = (q(a;) + q(ai41)/2. When the discontinuity crosses an adjacent cell center, 
the cell is updated according to the translated FWT curve fit give by Eq. 9. (Again, this update 


algorithm is simpler than that described earlier;'!? however, the results produce nearly identical 
results in the steady state.) 
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Figure 18. Solution to Burgers’ equation as a function of y and number of elements n with p = 4 (16 
cells per element) and m = 4 (fourth-order FWT fit) including shock detection. 


Results for the Burgers’ equation simulations are shown in Fig. 18. Symbols in Figs. 18(a- 
c) indicate cell center locations. Three diffusivity coefficients covering three orders of magnitude 
(v = 0.100, 0.010, and 0.001) are considered. When Az/v < O(1), the simulated shock is well 
resolved and when Az/v > O(1), the simulated shock is unresolved. All of the vy = 0.100 simulations 
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are well resolved. For example, Fig 18a, shows a well resolved shock with 8 elements and 16 cells 
per element. The LZ; norm for v = 0.100 in Fig. 18e shows a monotone decrease with slope less than 
-4, indicating fourth-order accuracy. The v = 0.010 simulations for n < 16 have Axv/yv = O(1). In 
these cases, the exact solution has only a few cells resolving the jump and the simulated solution 
either detects a shock with a jump across a single cell (as in Fig. 18b) or it places a discontinuity 
free FWT fit across the transition that tends to be sharper than the exact solution. The detection 
threshold (h*)!° for a discontinuity is set to 0.5 for these simulations; larger values tend to favor a 
smoother transition but may allow some small overshoot in the simulated profile. As the number of 
elements n increases, the value of Az/v becomes less than 1 and fourth-order accuracy is evident 
in Fig. 18e. The vy = 0.001 simulations for n < 32 have Ax/v < O(1). In this case, both the exact 
solution and the simulated solution show a jump across a single cell as shown in Fig. 18c. When 
n > 32, the ratio Ax/v becomes less than order 1 and the transition is well resolved, as shown 
in Fig. 18d. The Z; norm shows fourth-order accuracy in Fig. 18e for this range with n > 32. 
There is no clearly defined order of accuracy for smaller n; very small errors in the detected jump 
magnitude tend to keep the error norm of order 10~?. Numerical experiments indicate divergence 
for cases where the FWT polynomial has degree m > 4 and Az/v >= O(1). As with the advection 
equation, it appears that the extra degree of freedom in the FWT fit enables waviness across the 
detected discontinuity that amplifies with increasing time. 


C. Quasi-One-Dimensional, UnsteadyFlow 


The quasi-one-dimensional nozzle and shock tube problems are defined from the model system (Eq. 
17) with q = [pA, puA, pE A)’, f = [pwA, pu? A + pA, pull Al", h= vAzS4, and s = (0, =p24,0)" 
where primitive, dimensionless variables are density (p), velocity (u), pressure (p), and internal 
energy (e). The total energy (EF) is defined E = e + u?/2 and total enthalpy (H) is defined 
H = E+p/p. The equation of state for a perfect gas is p = Gpe where the ratio of specific heats 


for air as a perfect gas is y = 1.4 and 6 = y— 1. The speed of sound is c = \/yp/p. The matrices 
in Eqs. 31 - 33 are defined Aj41/2 = [u—¢,u,u+ rays and 


1 1 1 , sw +%é -é-Bi B 
Ri=| @-@€ wt w+é Ri = 38 22 -— pi? 2% —-28 
H-ié @/2 H+4é |, qoute CBG BO laa, 


A superscript * indicates a dimensional quantity. Pressure and enthalpy are nondimensionalized 
by the left state initial conditions; consequently p,es = 1 and Hef = 1. Density is nondimension- 
alized by p* = piep/Hrep 80 that preg = y/(y — 1). Velocity is nondimensionalized by ,/H7. ,. 

The detection threshold’® (h*) for a discontinuity is set to 0.3 for these Euler equation simula- 
tions. The acceptance threshold for a shock (see Fig. 13) is (pnigh/Plow)thr = 1.05 and for a contact 
discontinuity is (~pign/Plow)thr = 1.15. If these parameters are set too high, then the algorithm 
will fail to detect some weaker but real discontinuous flow features. If these parameters are set too 
low, then false positives for the presence of a discontinuity may be produced that tend to degrade 
accuracy and/or stability. 


1. Nozzle 


The nozzle area distribution A(x) used herein is defined 


A(z) = Athroat +3 (1 — cos(r2)) -l<a<l (45) 
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A nozzle schematic with boundary conditions is shown in Fig. 19. Two inflow boundary conditions 
at © = %q = —1 match total enthalpy and isentropic flow from the plenum. A single outflow 
boundary condition at x = x, = 1 matches pressure at the exit plane. 


A 


exit 


A 


throat 


p/p’ = Prien! Ppien 
-1 -0.5 0 05 x 1 


Figure 19. Schematic of nozzle geometry and boundary conditions. 


Analytic solutions to nozzle flow can be defined as a function of area ratio.2° The nozzle has 
Athroat = 1 and Aggie = 7. If 1 > perit > 0.9951909, then the flow throughout this nozzle is 
subsonic. If 0.9951909 > pexit > 0.180148, then the flow is sonic at the throat and a shock stands 
between the throat and the exit. The initial condition for all cases corresponds to a diaphragm 
at the nozzle throat separating plenum conditions from exit conditions with zero velocity and 
Dent orca)” = Dexie| Drek 

The test problem, with results at different times plotted in Fig. 20, sets p/prep = 0.2, y = 1.4 
and produces a standing shock near the exit with Mach number, M = 3.438 and a pressure ratio 
across the shock equal to 13.6 in the steady-state limit. The lines in this figure indicate the 
analytic steady-state limit solution. The symbols indicate cell center values from the simulation. 
The initial profile (Fig. 20 a) spawns a left moving expansion and right moving shock shortly after 
the diaphragm bursts (Fig. 20 b). A clearly discernible contact surface is not evident in the density 
distribution for this condition. At t = 0.646 the initial shock is near the exit plane and a steepening 
of the flow profile is evident in Fig. 20 c. In the next frame (Fig. 20 d), the steepening profiles 
transform to a right moving shock. The jump is recognized and preserved across a single cell by the 
FWT algorithm. By t = 2.66 the initial shock has exited, an expansion moves back up the nozzle 
and the shock “born” upstream in Fig. 20 d has moved far downstream, just past its steady-state 
position in Fig. 20 e. As time progresses, the shock moves back upstream and settles at the correct 
steady-state position in Fig. 20 f. 

The simulation in Fig. 20 used 16 elements across the domain (n = 16) with 16 cells per element 
(p = 4) and an FWT derived polynomial fit of degree 4 (m = 4). The steady-state solution for 
simulations n = 8,16,32 had associated error norms L; = 2.06 1074, 8.31 107°, 3.22 10-%°, 
respectively. This reduction in error norm indicates an order of accuracy equal to 1.36. The 
current algorithm makes no attempt to adjust the difference stencil based on the exact location of 
the discontinuity between the cells and so the net error is substantially weighted by this first-order 
treatment. Also, it was observed that a steady-state solution could be obtained without limiting the 
application of Eqs. 25 and 38 but the standing shock formed several cells upstream of the correct 
location. Consequently, this simulation was executed using Eq. 26 and a first-order extrapolation 
(L(a, b,c, d) = 6) in place of Eq. 38 whenever the stencil spanned a discontinuity. 
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(d) Shock formation near « = 0.3, 
t = 0.664 


(e) Standing shock at position closest 
to exit, t = 2.66 


(f) Steady-state solution, t = 50 


Figure 20. Evolution of the nozzle flowfield. Diaphragm at throat separates plenum conditions 
(p/Prep = 1, H/Hyes = 1) from exit condition (p/p,er = 0.2). The exit density is initialized assuming an 
isentropic expansion to the exit pressure from the plenum conditions. Simulation parameters in this 
case are: At/trep = 0.1Ar, n= 16, p= 4, and m= 4. 


2. Shock Tube (Sod Problem) 


The shock tube has closed ends and a diaphragm at the midpoint that separates the left state 
(p/Pref =1, p/ Pref = 6) from the right state (p/p;er = 0.001, p/ prep = 1). An analytic solution to 
this problem following the diaphragm removal is easily derived from an exact Riemann solver up 
to the first reflection of the shock off the right closed end. This test case (frequently cited as the 
Sod problem?°) produces a challenging contact discontinuity that follows the shock to the right. 
The simulation is executed with identical numerical parameters used in the previous nozzle case. 
The first three frames of Fig. 21 show the profiles at t/t;es = 0, 0.39, and 0.781, respectively. The 
lines show the time-accurate analytic solution assuming no reflections off the end walls. The symbols 
show the simulation results at the cell centers. Of particular interest is the density distribution 
shown by the red circles and dashed red line. The simulation preserves moving discontinuities at 
both the shock front and the contact discontinuity. Shock capturing techniques typically smear the 
shock over two to three cells and the contact discontinuity over many cells. The shock and the 
contact surface discontinuity in this simulation are generally within one cell of the exact solution. 
The error norm to the first reflection is shown in Fig. 22. The jaggedness of this plot is a consequence 
of the “quantized” movement of the discontinuities. The simulated discontinuity abruptly jumps 
from one cell interface to the next according to the logic in Fig. 14 whereas the analytic solution 
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moves continuously across the domain. 
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(i) t/trep = 14.454 


Flow in closed tube after diaphragm at midpoint is removed. Diaphragm separates 


conditions at right (p/pres = 1, p/pres = 6) from conditions at left (p/prer = 0.001, p/prer = 1). Simulation 


parameters in this case are: At/trep =0.1Axz, n = 16, p= 4, and m= 4. 


Figures 21 d and e show the reflected shock just after reflection and then shortly after it passes 
through the surface contact discontinuity. Note the contact discontinuity in the red circles at the 
extreme right of Fig. 21 e. The shock in this figure near x = 0 is moving to the left. 

The last four frames of Fig. 21 (fi) show a shock, initially at x = 0.6, moving to the right 
toward a contact discontinuity at « = 0.85. In frame (g), the oncoming shock has passed through 
the contact discontinuity at z = 0.9 and has just contacted the right wall. A weak reflected shock 
off the contact surface is transmitting to the left at « = 0.75. In frame (h), the strong shock has 
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Figure 22. Error norm as a function of time for the shock tube problem, up to the first reflection. 


reflected off the right wall and is overtaking the weak shock moving to the left near x = 0.35. In 
frame (i), the strong shock has overtaken the weaker shock and a single wave continues movement 
to the left at « = 0. 

The purpose of this test is to evaluate how well the FWT algorithm tracks shocks and contact 
discontinuities in a highly unsteady environment that includes colliding features. The use of analytic 
relations to define and move discontinuities enables preservation of even weak jumps for long periods 
of time. It is observed that the target norm for accepting the FWT fit (Eq. 27) is not always 
achieved when multiple discontinuities occur in the same element. This condition is observed in 
the time shortly after the diaphragm is removed and in the reflection of a strong shock off the right 
wall and in the vicinity of the near stationary contact discontinuity. While the majority of elements 
are well resolved by the FWT, it appears that the algorithm for detecting two discontinuities in 
the same element requires more refinement. 

In the current algorithm, conservation may be compromised every time a floater passes over a 
cell boundary. While the baseline formulation in smooth regions is conservative, there is currently 
no check that conservation is maintained when the floating boundary passes across a cell wall. 
While the exact Riemann solution is conservative, its implementation in discretized space requires 
additional constraints when the discontinuity passes over a cell center. These constraints are 
currently being investigated. A simple check of mass conservation in the current example should 
show that the integral of density with respect to distance from the contact discontinuity to the right 
wall should be a constant equal to 1.0. The integral of density over the entire domain, approximated 
by Sain piAx, should be a constant equal to 7. This integral over the entire domain varied between 
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6.78 and 7.07 (-3.1% to 1% error) over the time span of Fig. 21. (The same metric applied to the 
analytic solution up to the first reflection shows an error varying between -0.7% and 0.3%.) The 
integral of density from the contact surface to the right wall in Fig. 21(i) is 0.9675 (-3.25% error). 


VIII. Concluding Remarks 


The principal contribution of this paper is the evaluation of a Fast Walsh Transform (FWT) 
curve fit to enable feature detection in the solution of nonlinear differential equations. The FWT fit 
is enabled through a new theorem, Walsh Root Function Annihilation of Low-Order Polynomials. 
That is, given a uniform distribution of 2? cells on a one-dimensional element, it is proven that the 
inner product of the Walsh Root function for group p with every polynomial of degree < (p— 1) 
across the element is identically zero. It is also proven that the magnitude and location of a 
discontinuous jump, as represented by a Heaviside function, are explicitly identified by its Fast 
Walsh Transform (FWT) coefficients. These two proofs enable an algorithm that quickly provides 
a Weighted Least Squares fit to distributions across the element that include a discontinuity. It 
is shown that flux reconstruction relative to the FWT fit in partial differential equations provides 
improved accuracy and eliminates the need for flux limiting in the vicinity of a discontinuity. The 
detection of a discontinuity further enables analytic relations to locally describe its evolution and 
provide increased accuracy. 

Examples of inviscid flux formulation with respect to the FWT curve fits are provided for 
time-accurate advection, Burgers’ equation, and quasi-one-dimensional Euler equations. An RK-4 
method is used to advance all solutions in an explicit, time-accurate formulation. Fourth-order 
accuracy is confirmed for the advection of a Gaussian profile and for the steady-state solution of 
Burgers equation with a cell Reynolds number Az/v < 1. Discontinuities in a test profile translate 
across the domain twenty times without dissipation (jump across a single cell) and without any 
overshoot, maintaining second-order convergence. Discontinuities in the nonlinear sample problems 
are detected within a single cell, without ringing using the FWT fit formulation. The steady state 
limit for a standing shock in a de Laval nozzle is accurately recovered. The algorithm is able 
to detect and preserve shocks and surface contact discontinuities through multiple reflections and 
collisions within a closed, shock tube problem. 

Some thirty years ago, the CFD community predominantly adopted upwind-based, shock- 
capturing methods over various implementations of shock-fitting methods, utilizing unsteady Rank- 
ine Hugoniot relations to compute analytic jump conditions and shock velocities. Shock-capturing 
algorithms make minor adjustments to a stencil based on local gradients to enable a stable dis- 
cretization across the shock (limiting) but do not use any analytic information based on the Rankine 
Hugoniot equations. Shock-fitting algorithms must carry the infrastructure to define interior shock 
coordinates and make significant adjustments to the stencil at these interior boundaries, but they 
benefit from a more accurate, analytic representation of shocks moving through the domain. The 
vision for the FWT fit is to work like a shock-capturing method with respect to detection but work 
like a shock-fitting method with respect to use of more accurate analytic relations. The enabling 
innovation for feature detection is the ability to quantitatively define the magnitude and location 
of a discontinuity within some 2? cells in an element within the domain using the FWT curve fit. 
The cost here is a larger stencil of cells required than for either shock capturing or shock fitting. 
The expectation is that the state of the art is moving toward use of higher-order elements so that 
the requisite stencil for detection is already available. 
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